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ABSTRACT 



CN ■ The multiple images of lensed quasars provide evidence on the mass distribution 

of the lensing galaxy. The lensing invariants are constructed from the positions 



of the images, their parities and their fluxes. They depend only on the structure 



00 
O 

C\| ■ of the lensing potential. The simplest is the magnification invariant, which is 

2 ' the sum of the signed magnifications of the images. Higher order configuration 

O ■ invariants are the sums of products of the signed magnifications with positive or 

^ ■ negative powers of the position coordinates of the images. 

I ■ We consider the case of the four and five image systems produced by elliptical 

2 ■ power-law galaxies with tjj oc (a;^ + This paper provides simple contour 

^ ■ integrals for evaluating all their lensing invariants. For practical evaluation, this 

offers considerable advantages over the algebraic methods used previously. The 
magnification invariant is exactly B = 2/(2-/9) for the special cases jS = 0,1 
and 4/3; for other values of f3, this remains an approximation, but an excellent 
one at small source offset. Similarly, the sums of the first and second powers 
of the image positions (or their reciprocals), when weighted with the signed 
magnifications, are just proportional to the same powers of the source offset, with 
a constant of proportionality B. To illustrate the power of the contour integral 
method, we calculate full expansions in the source offset for all lensing invariants 
in the presence of arbitrary external shear. As an example, we use the elliptical 
power-law galaxies to fit to the data on the four images of the Einstein Cross 
(G2237+030). The lensing invariants play a role by reducing the dimensionality 
of the parameter space in which the minimisation proceeds with consequent 
gains in accuracy and speed. 



Subject headings: gravitational lensing - galaxies: structure - quasars individual: 
G2237+030 
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1. 



Introduction 



There are now some sixty or so multiply imaged quasars known. Most of these exhibit double 
or quadruple images, although some triples and rings are also known (see Warren et al. 1999, 
Ratnatunga, Griffiths, Ostrander 1999, Wisotzki et al. 1999 and Tonry & Kochanek 1999 
for some recent examples, as well as Pospieszalka et al. 1999 for details of the gravitational 
lensing database which maintains a list of candidates) . This paper is mainly concerned with 
the sixteen or so quadruple systems. In propitious circumstances, the observables of such 
systems are the location of the lens, the positions of the four images, their parities and 
their fluxes. Lensing invariants are constructed directly from the lensing observables and 
remain invariant against changes in the source or lens configuration (provided a caustic is 
not crossed). A classical problem in the theory of gravitational lensing is the construction of 
a lens model that reproduces the observed properties of the images. One of the advantages 
of lensing invariants is that they provide simple tests by which we can determine quickly 
whether a given set of images can be produced by a particular lens. 

Witt & Mao's (1995) paper provides the ffi^st example of a lensing invariant. They 
considered a binary lens of two point masses, which has either three or five images according 
to whether the source lies outside or inside the caustic (e.g., Schneider, Ehlers & Falco 1992, 
section 8.3). Their main result is that the sum of the signed magnifications of the five images 
is unity when the source lies within the caustic, namely 



Here, jj,i is the absolute magnification of the image, while Pi is the parity. This result holds 
good irrespective of the position of the source, provided it remains within the central caustic. 
Subsequently, Rhie (1997) showed that the sum of the signed magnifications of the images of 
the N point mass lens is also equal to unity, provided the source location yields the maximum 
number of + 1 images. These are remarkably simple results, bearing in mind that the 
lens equation that provides the locations of the images must be solved numerically. Further 
examples of lensing invariants were found by Dalai (1998) and Dalai & Rabin (2000). 

Witt & Mao's (2000) paper provides a study of lensing potentials stratified on simi- 
lar concentric ellipses and falling off like a power-law with index /3 suitable for modelling 
quadruple lenses. They showed that for the point mass (/? = 0) and isothermal (/3 = 1) 
cases, there is an invariant 



They introduced higher-order moments and showed that they too had simple results in the 
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point mass and isothermal cases 

4 4 

^ HiPiXi = B^, ^ HiPiVi = Br]. (3) 
i=i i=i 

Here, (xj, yj) are the coordinates of the ith image, whereas (^, 77) is the position of the source. 
Witt & Mao (2000) argued on the basis of numerical experimentation that (2) and (3) remain 
good approximations for all other values of (3. They suggested the use of interpolation to 
extend the result to all the power-law potentials. Witt & Mao's (2000) paper is especially 
important to us here, as our results extend the work that they began. 

In §2, we exploit contour integral methods to provide a new way of finding the invariants 
of a gravitational lens system. The idea is very simple. The images are the roots of the lens 
equation. So, by a judicious choice, the observablcs at each image can be made to equal 
the residue of a complex integrand. The lensing invariants, which are sums over the images, 
then become equal to a contour integral by Cauchy's theorem. Although our method is more 
general, this paper concentrates on applications to the elliptic power-law potentials. In §3 
we evaluate invariants for the case in which S = 2/(2 — /3) is a positive integer. Our complex 
integrand is then single- valued. The simplest cases are the point mass (/3 = 0, S = 1) and 
the isothermal (/3 = 1, 5 = 2) lenses noted by Witt k Mao (2000). The invariants for S = 3 
include the contributions from a weak fifth image, while those for integer B > ?> include small 
contributions from spurious images, which are non-physical solutions of the lens equation. 

In §4, we progress to the case of non-integer B when the contributions from a branch 
cut must be added to the earher results. Remarkably this gives us exact expressions for 
the lensing invariants for the four bright images for all B, because the added contributions 
cancel out the effects of the fifth and the spurious images. Thus, our contour integral method 
provides not only a simpler but also a more powerful way to determine the lensing invariants 
than the earlier algebraic methods. Astrophysical applications of the invariants are presented 
in §5. They are useful because they can short-cut the modelling process, telling us quickly 
whether an assumed lensing potential (in this case the elliptic power-law potentials) can 
reproduce the data. We pay particular attention to the Einstein Cross gravitational lens 
system (G2237-I-0305) as our example. Finally, §6 discusses our conclusions and presents 
plans for future work. 
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2. The Contour Integral Method 

2.1. The Lens Equation 

At outset, we shall assume that the lensing potential is stratified on similar concentric 
ellipses 

^{x^ + y^q-y/^ if0</?<2, 

(4) 

4log(x2 + yV') if/9 = 0. 



This family of models was introduced into gravitational lensing by Blandford & Kochanek 
(1987) and subsequently studied by others (e.g., Kassiola & Kovner 1993; Witt 1996; Witt 
& Mao 1997, 2000; Evans & Wilkinson 1998). Here, q is the axis ratio of the projected 
equipotentials and is chosen to satisfy < g < 1 without loss of generality. The parameter 
(3 controls the radial profile of the potential while A measures the depth of the well. The 
point mass (Schwarzschild lens) has j3 = and q = 1, while the isothermal sphere has /3 = 1 
and q = I. The gravitational convergence (or projected density) is 

_ A [l + q^{p-l)]x^ + [l + q-^((3-l)]y^ 

This is positive definite provided (3 > 1 — q^. Three-dimensional analogues of these models 
are already famihar in galactic dynamics as the power- law models (Evans 1993, 1994). They 
are the only fiattened and reasonably realistic galaxy models known, for which the self- 
gravitation equations can be solved to find simple two-integral distribution functions. For 
example, they have been used in modelling the nearby elliptical galaxy M32 (van der Marel 
et al. 1994), the inner parts of the Galactic bulge (Evans & de Zeeuw 1994), as well as the 
dark halo of our own Galaxy in the interpretation of the microlensing results (Alcock et al. 
1997). 

For the geometrically thin lens with projected potential (4), the paths of photons are 
given by (e.g., Schneider, Ehlers & Falco 1992, section 5.1) 

Ax MlQ~^ 
e = a;-h7iX + 72l/- ^^,^^,^_,^,_^/, , = ^ + ^'"^ " " (^^T?Pj^' ^ ^ 

where (^, rf) are Cartesian coordinates of the source. Here, 71 and 72 allow for a constant 
external shear in an arbitrary direction. The external shear may be produced by the gravity 
field of a cluster or by the effects of nearby galaxies and usually hes between < I7I < 0.3 
(Keeton, Kochanek & Seljak 1997). We shall always assume that the lens is not circularly 
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symmetric (that is, either 1 or 71 or 72 7^ 0). The lensing properties of circularly 
symmetric lenses are qualitatively different from those of more realistic flattened lenses. For 
example, the tangential caustic degenerates to a point in the case of circular symmetry (e.g., 
Schneider et al. 1992, section 8.1). So, circularly symmetric lenses can give misleading 
results and they are not widely used in modelling. 

A formulation in terms of complex numbers has been used before to ease calculations 
in lensing theory (e.g., Bourassa, Kantowski & Norton 1973, Bourassa & Kantowski 1975, 
Witt 1990). However, we shall find it helpful to use somewhat different definitions, namely 

C = C + m, z^x + iy/q. (7) 

The lens equation (6) becomes 



Az 

C = |[1 + + 7i(l - q')]z + i[l - + 7i(l + q') + 21^72]^ - 7— OT- (§) 

[zz) ^' 



It is convenient to introduce a variable t such that: 



_ A 2 
^~(ii)V^' ^~2^' 

where the bar represents complex conjugation, and henceforth to set the constant A equal to 
unity. The dependence of our results on A can be recovered by dividing all powers of z and 
( and their conjugates by the same powers of A^^'^. Then, the lens equation can be written 
in matrix form as 

with 

P^l[l + q^ + 7i(l - _t^p^-t, Q^l[l-q^ + ^,(1 + g^) + 2iqj,], (10) 
from which it follows, by the standard algorithms of matrix inversion, that 



Z J P2_ |g|2 y_Q P J \ C 

By forming = l/{zz), we see that the solutions of the lens equation also satisfy the 
equation 



K{t;C,0 = t''{PC-QO{PC-QO- 



P'-\Q\' 



0. (12) 



Let us designate K{t, C, C) =0 as the imaging equation. Its real and positive solutions 
provide the possible image positions t (or equivalently z), given the source location C and 
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its conjugate (. The imaging equation must have at least as many solutions as the original 
lens equation (6). It may have more; that is, there may be solutions of (12) that do not 
correspond to true images and which we will call spurious roots. 

For example, when B is an integer N < 2, then the imaging equation is a quartic 
with four roots each of which corresponds to one of the images of a quadruple lens. When 
B = N > 2, then the imaging equation is a polynomial of degree N + 2. In the case B = 3, 
the roots all correspond to the five images of a quintuple lens. However, when B > 3, there 
are always roots of the imaging equation that are spurious, as there are never more than five 
images. This is shown in Appendix A. 



2.2. The Magnification Inveiriant 



Mathematically speaking, the lens equation defines a map from the lens plane to the 
source plane (see Schneider et al. 1992, chapter 5). The Jacobian matrix of this mapping is 
denoted by J, where 



J 









/ m 


z W 








) 


\ m 


z ^ 


z ' 



J- 



/ dz 




dz 




{ K 


c 


dC 




dz 




dz 


;) 


\ K 


c 







With our choice of source and image variables, the reciprocal of the determinant of J corre- 
sponds physically to the signed magnification of an image via 

•,2 

fJ^iPi, (13) 



detJ 



where //j is the absolute value of the magnification and pi is the parity of the image located 
at {xi,yi). For certain positions, detJ may vanish and the magnification is infinite. These 
are the critical points and lines. The caustics are the images of the critical points and curves 
under the lens mapping (6). 



By differentiating the expression for t, we obtain: 



B dt 



t^+^ dC 
B dt 



t^+^ dC 
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Recalling the properties of the Jacobian matrix J and its inverse J , it follows that 
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(14) 
(15) 

(16) 
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dz 



_H — 3 
C dz 



^l+B 



-z. 



;i7) 



By multiplying (16) by and (17) by and subtracting, we construct the determinant of 
the Jacobian matrix and so establish 



c 



dz 



B 



dz 



(18) 



However, equation (12) provides us with the requirement that K{t; C) C) — and so, using a 
standard formula in the theory of partial differentiation, we have 



1 



B 



detJ 



'dC 
'Uz 



dK 




dC 




dK 




dt 





(19) 



Now comes a remarkable simplification. For the ellipsoidally stratified potentials, we find by 
direct evaluation that 



dC 
dz 



dC 
dz 



[P^ + \Q\^)C - 2PQC 



A similar factor occurs on partially differentiating K{t; C, C)j other words 

dK 



(20) 



dC 

So, equation (19) simplifies to 



{P' + \Q\')C-2PQC 



B P^-\Q\' 



detJ 



t dK 



(21) 



(22) 



Let us recall that the images correspond to simple zeros of K{t;(X) and the residue of 
1/K at a simple zero is just l/Kf. This gives us a contour integral representation of the 
magnification invariant, which is defined as the sum of the signed magnifications of the 
images 

q^B j;dtP^-\Q\^ 



images 



27ri Tc t K{t; C, C) ' 



(23) 



where C is a contour in the complex lens plane enclosing only the simple poles corresponding 
to visible images. It is an invariant, as it is a pure number which depends only on the 
structure of the lensing potential or, equivalently, the mass distribution. 
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2.3. The Configuration Invariants 

The higher order invariants are sums of powers of the coordinates of the images weighted 
with their signed magnifications. We shall call them the configuration invariants. From the 
lens equation (6), we deduce that 

The higher order configuration invariants are the sums over the images of products of the 
signed magnifications with the position coordinates. The mth moment with respect to x and 
nth moment with respect to y is given by the contour integral 

E ^^^P^xTy: = / f C, OrW, C, C)]". (26) 

images \ iSiS) 

Witt & Mao (2000) introduced some higher order moments of x and y for the cases when 
the integers m and n are positive and they showed that the results for the lensing potentials 
(4) were very simple for the point mass (P — 0) and isothermal (/? = 1) cases. 

It is also possible to derive moments with respect to the complex image coordinates. 
The mth moment with respect to z and the nth moment with respect to z is given by the 
contour integral 

ii^es " 27ri % t Kit; (, C)[P^ - IQIT^^-' ^ ^ 

For positive moments, the results are of course just a recasting of (26), but for reciprocal 
moments, the results are different. In both cases (26) and (27), the contour C is chosen 
to enclose the poles corresponding to the visible images, but to exclude other poles in the 
complex i-plane such as those at the zeros of — when m + n > 1. 



3. Invciriants for integer B 

To begin with, let us study elliptic power-law potentials for which the exponent B — 
2/(2-/3) is an integer. This may strike the reader as narrow in its scope, but it is not 
really so. As we will see in §4, results which are exactly true for these cases are also good 
approximations for the four bright images of all the power-law galaxies if the source is 
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sufRciently close to the galactic center. The simplest cases are B — 1, or (3 — 0, when the 
potential approximates that of a point mass plus shear, and S = 2, or /3 = 1, when the 
potential is isothermal-like and represents a galaxy or cluster with a fiat rotation curve. In 
both these cases, the imaging equation (12) is a quartic and has four roots. When i? = 3 or 
(3 = 4/3, so that the convergence or projected mass density falls off like r^^/^, the imaging 
equation (12) is a quintic and has five roots and there is a weak fifth image. The cases 
B — 1,2 and 3 are special because there are then no spurious roots of the imaging equation 
when the source lies within the central caustic. There are spurious roots of the imaging 
equation for S > 3, which contribute to the invariants which we calculate. However their 
contributions are small and many of them cancel, with the consequence that the results 
remain useful even when there are spurious roots. Cases of integer B avoid the need for a 
branch cut to define the term in the imaging equation (12) as a single- valued function in 
the complex plane. 



3.1. The Magnification Invciriant 



As the images are given by the zeros of the imaging equation, we have 

q^B fdtP^- |gp 
2^ JcT K{t;C,C) 



q^B rdtP^-\Q\^ 



images 



The integrand has simple poles at the zeros oi K{t;(X)^ ^ well as at the origin. The contour 
C is traversed in an anti-clockwise direction and encloses the poles corresponding to images. 

Suppose the integrand is taken over a large circle at infinity Coo- As i — > oo, the 
integrand falls off at least as fast as 0{t~^), so we have; 

. ^^^^0. (29) 

Distorting the contour so it encloses the zeros of K{t; C, C) ^^'^ the pole at the origin, we 
have 

V- q^B r dtP^-\Q\^ , , 

> y^iVi + 7^ <P , = 0, 30 

^ 2;ri Jc t K(t;CX) 

images V ) -b ? -b / 

where is a small circle of radius e enclosing the origin in the anticlockwise direction. As 
t 0, then K{t; C,C) -{P^-\Q\'^f, where Pq = P{t = 0) = |[1 + + 7i(l - g^)] so that 

^ _ Bq' _ B 

images " ' ' ' ^ ' ^ 
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This remarkably simple result was first obtained by Witt & Mao (2000), who deduced it for 
the B = 1 and 2 cases and used interpolation to extend it as an approximation to all B. 
In fact, the result (31) is exact for S = 3 as well, provided the contribution from the fifth 
de-magnified image is included. 



3.2. The Configuration Invariants: First and Second Moments 

The first moments can be obtained in a similar way. When m = 1, n = the integrand 
(27) still has poles only at the zeros of the imaging equation and the origin. Evaluating the 
residue at the origin, we find 

^ _ ^[PqC-QC] ...^ 

images ^ ^ "- 

and hence 

images ^ '-^ images ^ ' ' 

Witt & Mao (2000) give these results in the case of on-axis shear only (72 = 0), while Dalai 
& Rabin (2000) give them for the case B = 2 only. Algorithms using the real analysis tend 
to become cumbersome with off-axis shear. It is part of the power of the contour integral 
method that results can be obtained easily for general shear. When there is no off-axis shear, 
the X (or y) moments of the signed magnifications of the images are related only to the ^ 
(or rj) offset of the lens. 

The second moments are given by m + n = 2. The integrand (27) now has additional 
simple poles at the two real points ti and t2, where 

ti^Po+\Q\, t2^Po-\Q\, (34) 

where P^-\Q\'^ = 0. Provided that 72+7! < l, then Pq > and P^-\Q\^ = g^(l-7?-7|) > 
0, and both ti and ^2 are positive with ti > t2 and ^2 < 1- The residues at these locations 
give rise to constant terms independent of the source position in the configuration invariants. 
To show this, let us consider the evaluation of one of the second moments (m = 1, n = 1) in 
a little detail. At large radii, the integrand falls off at least as fast as 0{fr^). If evaluated 
about a circle at infinity Cqo, the contour integral (27) vanishes. Distorting the contour so it 
encloses the zeros of K{t;(X) and the poles at t = 0, ti and t2, we have 

- ^<fB r dt \P(.-QZ\[PC,-QC\ „ 

,£''*-+2^i^..,.,T[P^-|«Pli^(t;C.C)^''- 
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Here, Ce, Ct^ and Ct2 are small circles of radius e enclosing the simple poles at t — 0,ti and 
t2 respectively. By evaluating the residues, we find; 

V a-p-z-z- - i!^ / J L\ + B(PoC-QOiPoC-QO 

- 2\Q\ U-^^ tf-^^ S + ,^(1 - 7? - 7i)^ ^ ^ 

After evaluating the moment of zf in a similar manner, we obtain the real second moments: 



images 



^ 2 _ q'B [ |Q|+Re(g) |g|-Re(Q) ) ^ 



images 



■.hs ~ 4|g| + (i-7f-7i)^ ■ ^ ^ 

These results simplify considerably when the shear acts only on axis (72 = 0), in which case 
\Q\ = Ke{Q). ti = 1 + 7i, and t2 = q'^{l — 71). The (or xy or y^) moments of the signed 
magnifications of the images are related only to the (or ^rj or 77^) offset of the lens. The 
results (37) - (39) are given by Witt & Mao (2000) for the case of vanishing off-axis shear. 
For completeness, the third moments are derived in Appendix B. 



3.3. The Configuration Invciriants: Reciprocal Moments 

The reciprocal moments can also be found. The first order reciprocal moments are given 
bym = — l,n = Oorm = 0,n = — 1. In addition to the poles at the images and the origin, 
the integrand also has simple poles at either ^3 or given by 

ts = ?^(l-Ti-^), «. = 1+Ti-^. (40) 

The real first order reciprocal moments are 

The constants Ci and C2 are 

q^B^I _ 57| 



mX-ll] [(l-7i)C-72r?][-7l + e^-^g'^((l-7i)e-72^)^]' 

Bll 

r]h[r]^t^ - ll\ [(l+7i)^-72e][-7i + ^'-^((l+7i)^-720^] ' 
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Sending ^ ^ r], r] ^ ^, ^ —71, — > 1, then we see that Ci — > C2. This follows as Ci 
is the residue of the integrand at ta, and C2 at t4. Note that Ci and C2 both vanish when 
the shear is on-axis only (72 = 0). In this instance, it is again true that the reciprocal x (or 
y) moments of the signed magnifications depend only on the reciprocal ^ (or rj) offset of the 



source. 



The complex reciprocal moments have poles of order — m or —n dX t = = Pq — QC/C 
and/or its conjugate t = = Pq — QC/C- For m — —1 and n — 0, the residues at t and 
t — exactly cancel so that 



0. 



(44) 



images 



This implies the vanishing of the following two real moments 



E 



xf + yfq- 



E 



xf + yfq- 



0. 



(45) 



images ' " images 

A complication that arises with sufficiently negative m + n is that there is an additional 
contribution from the integration over the large circle at infinity C^o to be included. This 
happens when —{m + n) equals or exceeds the larger of B and 2, and hence with second 
order and higher reciprocal moments. The second order reciprocal moments in complex form 
have the compact expressions 



images 



fJ-iPi 



0, 



2g^ 



images 



-2\2 



Bq^ 



ICP 



B <2, 
B^2, 

B>2, 



(46) 



and 



y2 — 



^^ipi 



images 



images 



{x\ + ylq-'^f 



2ixi?/i 



Bq^Qi 



2q^ 



Bq' 



_l 

ICP-i \Q? 
1_ _ Qt\ 



B<2, 
B^2, 
B>2. 



(47) 

There are also moments for which one of m and n is negative, and the other positive. The 
simplest is that for m = 1 and n — —1, which gives 

BPo 



E 



images 



IJ'iPiZi 
Zi 



0(1 - 1! - 7l) ' 



(48) 
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or, in real form, 



E 



liiVm - m ) BP„Re(Q) 



3.4. Spurious Roots 

Now let us confront the problem of the spurious roots of the imaging equation (12) that 
do not correspond to physical images. When -B < 2, there are four roots on the real positive 
i-axis corresponding to the four images of a quadruple lens. When B > 2, there are B-\-2 
roots in total. With \C,\ small, the B — 2 additional roots appear at large t (or at small |2;|). 
By balancing terms in the imaging equation, we see that these roots satisfy to leading 
order 

^^ = p^-p(|?^)+o(i)' (50) 

where n is an integer so that there are B — 2 distinct roots of this form. Only one of them 
{n — 0) lies on the real positive i-axis and corresponds to a physical image. In fact, it is 
the strongly de-magnified fifth image that occurs so close to the centre of the lens that it is 
often hard to detect. The other B— 3 roots (0 < n < -B— 3) are spurious and do not provide 
physical images. Figure 1 illustrates the behaviour of the roots in the complex t-plane as B 
increases. The first panel shows the case 5 = 3, with the five roots corresponding to the 
five images along the real positive t-axis. The second panel shows the case 5 = 4, when a 
spurious root appears at large negative real t. When 4 < 5 < 6, this negative real root splits 
into a pair with arguments exp ±27ri/(i? — 2). In the third panel, at S = 5, we have three 
roots equally spaced around a circle of large radius in the complex i-plane. In the fourth 
panel, the cycle of events begins again with a new spurious root appearing on the negative 
real axis when B — Q. 

Some small part of the invariants that we calculated earlier in this section are due to 
spurious roots. Their contributions can now be calculated and subtracted. We shall do this 
for the sum of the signal magnitudes. In fact, we shall subtract also the contribution from 
the weak fifth image to obtain an approximation for the sum of the signed magnitude of the 
four bright images. The contribution to the integral (28) from the residue at each root (50) 
is 

Residue(t = Q = exp + 0{\Cf^^-^^ (51) 
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Fig. 1. — The locations of the roots of the imaging equation in the complex i-plane for the 
cases B — 3,4,5 and 6. New spurious roots appear on the negative real axis whenever B is 
an even integer. As B increases, the spurious root bifurcates and the pair move off into the 
complex plane with arguments exp ±27ri/(i?— 2). For integer B, the roots lie at the vertices 
of a regular B — 2 sided polygon. 
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For odd values of S, it is convenient to set B — 3 + 2N. The spurious roots lie at angles 
exp(±27rin/(2A^+l)) and so 



4 images 
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Bq^ 



(l-7?-7l) {B-2) 
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ICI 



4/(S-2) 
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n=l ^ ' ^ 



+ 0(|Cr/(^^)). (52) 



For even integers B — 2+2N, the spurious roots he at angles exp{±n7ri/N) and so 



B 



Bq^ 



4 images 
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n=l ^ ^ 



ICI 



4/(S-2) 



+ exp - 



Arni'i 
1^ 



+ 0(|Cr/(^^)). (53) 



Except for the cases B — 3 and B — A, the terms in square brackets on the right-hand side 
of (52) and (53) sum to zero, and so there is no contribution to the magnification invariant 
at this order. This is not the case for S = 3, when the sum in (52) is empty, and we obtain: 



E 



(1 - il - iD 



-3g^icr+o(icn, 



(54) 



4 images 

or for = 4 when the sum in (53) is empty, and we obtain 

4 



4 images 



(1 - il - iD 



-vicr+o(icr). 



(55) 



In both cases, the corrections match those given by the more general theory which we develop 
in the next section. 



4. Invciriants for all B 

The general case in which B is no longer an integer requires that a branch cut be 
introduced to make single- valued the term in the imaging equation (11). This apparent 
complication is a blessing in disguise. It allows us to generate invariants which are valid for 
a continuous range of values of B. They are infinite series in the source coordinates, whose 
leading terms are just those that we derived in §3 — namely, eqns (31) - (33), (36) - (39), 
and (41). Importantly, they are also analytic functions of B, and hence are useful for all 
values of B in the physically relevant range l<S<oo (0</3<2). This section explains 
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why this is so. We first derive the magnification invariant for the case of general B < 2 
(§4.1) and then the configuration invariants (§4.2). In §4.3 we show how these invariants 
could also be derived more laboriously by perturbation theory. Readers happy to take the 
results on trust can turn directly to the applications (§5). 



4.1. A Recalculated Magnification Invariant 

Suppose that B is not an integer and that 1 < S < 2, so that there are no spurious 

roots and only the four bright images. We make single- valued and real on the positive real 
t-axis with a cut along the negative real axis. We again use the contour integral (28) for the 
magnification invariant. We enlarge the contour C as much as possible — the integrand has 
no singularities other than those on the cut — until it consists of the infinite loop Ccut shown 
in Figure 2 together with a large circle at infinity. The large circle again contributes nothing, 
because of the decay of the integrand for large t. Hence, the contour integral becomes 

g'B f dtP^-\Ql 

Let us write K — Ki + K2, where 

K^ = -(P2 - \Q\y = -{t-Uf{t-h)\ K, = t^(PC - QC){PC - QC). (57) 

We introduce the following two linear combinations of the complex source coordinates: 

They have the important properties that 

Ai/ + AP = 0, |A|2 + = ICI^ (59) 

and allow us to write 

PC-QC = X{t2-t)+u{t,-t), (P(-Q()(P(-Q() = \x\'{t2-tf+\u\'{t,-tf. (60) 

We expand the denominator of our integrand in an infinite series to get 

This is vahd for sufficiently small |C| because both \Ki\ and 1^21 are bounded below on the 
cut for finite t and \K2/Ki\ — > as \t\ — > 00. The cut is unnecessary for the first j — 
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term of the series (61). It has just a simple pole ed, t — and so can easily be evaluated. 
The result is exactly the same as that of §3.1 as the residue calculation is unchanged. In the 
other terms of the integrand, we set t — ■u;exp(±7ri) on the upper and lower sides of the cut 
to get 



4 images 



l-7i^-7l 



+ 



.7=1 



Binomial expansion of the integrand gives us integrals which are all of the form 



sin Bjn 



w- 



1<B<2, 



(62) 



(63) 



TT Jo {w + hy+\w + t2y^-'+'' 

for integers i in [0,2_7]. The integrals can all be expressed in terms of Gaussian hypergeo- 
metric functions 2-^1, using formula [3.197.1] of Gradshteyn & Ryzhik (1965) as 



_1 2^'+^ 



.jB+£-2j-l 



s=l 



''1 



-2Fi(£+l,ji?;2j + 2;l-|). (64) 



Appendix C gives closed formulae for some of the low order integrals Ij^e. Binomial expansion 
of the numerator of the integrand gives the following infinite series for the magnification 
invariant of the four bright images. We henceforth label this as M.{B) to emphasise that it 
is an analytic function of B; 



M{B) 



4 images 



B 
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l-7f-7l 
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(65) 



The jth term is a homogeneous polynomial of degree 2j in the source coordinates. Writing 
out the first correction term explicitly, we have 



MiB) 
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7l 



+ Bq^ [K^.e + Kll^rj + KW] + OdCll, 



(66) 



where 
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(67) 



and, using eqn (C2), wc find, 
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(68) 
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This gives explicit formulae for the first order corrections to the magnification invariant. 
Once again, the formulae simplify considerably if there is no off-axis shear and Q is real. 

Although the preceding derivation of the infinite series (65) for the magnification invari- 
ant is valid only for 1 < S < 2, the series itself is not so limited. As we show in Appendix 
D, it converges for arbitrarily large values of B for sufficiently small In fact, the series is 
the unique analytic continuation of the magnification invariant M.{B)^ which is an analytic 
function of i?, and hence can be used for B > 2. Here are two simple tests that give credence 
to this remarkable result, the justification of which is given in Appendix D. First, for B = 3, 
Ii,2 = h,o = and there is no 0(|Cp) term in the series (65). At the next order, we find 
from equation (64) that /2,4 = l2,2 = -^2,0 = ~1 and hence 

M{3) = . J - 3g2|Cr + 0{\C\% (69) 
1-71-72 

in agreement with the result (54) of §3.4. Second, let us set B — A. Then /i_2 = -^1,0 = — 1, 
and evaluating the j — 1 term of (65) and using (59) gives 

M{4) = ^ - 4q'\C\' + O(ICr). (70) 

7i 72 

This matches the approximate result (55) we obtained in §3.4 after subtracting the weak 
fifth and one spurious image. 



4.2. Recalculated Configuration Invariants 



We now obtain infinite series expansions for the configuration invariants (27), also 
summed over the four bright images, by manipulating that contour integral in a manner 
similar to that of the previous section. The integrand now generally has additional finite 
singularities at one or more of ti, t2, h, and ^5, and so their residues must be accounted for 
when they are crossed in the process of enlarging the contour and wrapping it around the 
cut. The result is: 

B[Poc-Qcr[Poc-Qcr 



4 images 



m 



Q 



2m+2n 



[1- 



q^B ^ residues of 



poles 



tK{fXX)[P^-\Q\^]-^+-' 



+ Sm,n{B), (71) 



where 



Sm,n{B) = Bq^J2 
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(72) 
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Re(t) 



Fig. 2. — The contour for eq. (56) is shown. The t-plane is cut along the negative real axis 
to make single-valued. The contour is completed by a large circle at infinity, along which 
the integrand is vanishingly small. 



-20- 



The binomial expansion of the integrand is less simple than in §4.1, as m 7^ n and a double 
summation is needed. The result is the following homogeneous polynomial of A, A, and P: 



00 2jim\n 



fe=min(£,mr|-j) 



m+ A / n+i 



k 



■k 



yl-4i ^in\j--k -n\i^-4 (12)) 



The jth term is now a homogeneous polynomial of degree 2j+m+n in the source coordinates. 
A more general integral with an extra subscript (which was zero previously) now appears as 
a coefficient. It is 



Ij,e,N{ti,t2) — 



sin Bjir r 
TT io 
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(74) 



/_-|\Ar+l 2j4-iV+l jB+E-2j-N-l / , 

= p^TiW n Wi5-)^^.F.(£+UB;2j>2+7V;l-^ 

Although the integral definition of Ij,e,N is valid only if i? < 2, the hypergeometric formula 
is valid for all 5 > 1. Explicit closed form expressions for all the hypergeometric functions 
are given in Appendix C. Using these results allows us to extend the formulae given in §3.2 
and §3.3 to configuration invariants. The larger m + n is, the longer they are, so we shall 
here quote only the first reciprocal moments, which are 
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(75) 



The first order corrections for the first moments are given in Appendix C. 



4.3. Perturbation Analysis 



Invariants, and corrections to them, may also be derived by using perturbation expan- 
sions to solve the imaging equation. Its four main roots tend in pairs to ti and t2 as the 
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source tends to the center of the lensing galaxy, and the expansions are in integer powers of 
C and One pair of images hes at 



1 \|2 _i_ 1+3B/2| 



and the other at 
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+ 0(|Cr), (77) 



where A and v are as defined in equation (58). The signed magnifications for the first pair 
of roots are, from equation (22), 
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(78) 



while those for the second pair can be found from them by applying the transformations 
ti ^ t2, \ ^ u, and \Q\ — > — \Q\. Summation then gives equation (66) again, as it should. 

Finding configuration invariants by perturbation expansion requires the additional ef- 
fort of expanding for the image positions z and z. Though the individual terms are of 
some intrinsic interest, perturbation analysis is an unnecessarily laborious way of computing 
invariants now that we have the contour integrals which compute them to all orders. 



5. Astrophysical Applications 

There is inevitable uncertainty in the distribution of mass in lenses, and all methods 
(whether parametrised fitting or non-parametric modelling) must make strong assumptions 
as to the mass distribution before they can make progress. The major application of lensing 
invariants is to short-cut the modelling process. In §5.1, we provide simple tests to deter- 
mine whether a given lens system is well modelled by an elliptic power-law potential with 
shear. Next, we use the lensing invariants to develop an algorithm for estimation of the 
model parameters and apply it to fake data in §5.2 and real data for the Einstein Cross 
(G2237-I-0305) in §5.3. The work in §4 has demonstrated that the simple equations for the 
moments we derived in §3 - namely, eqns (31) - (33), (36) - (39), and (41) - are excellent 
approximations for all B. Here, we shall assume that these results are exact. 
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Fig. 3. — The recovery of the test ratios (81) in the cases of 5% and 15% errors on the 
photometry. The lens equation is solved to find the true image positions and fiuxes, then 
photometry errors are added and the test ratios computed. The histograms show the distri- 
butions obtained with 400 realisations. 
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Fig. 4. — The recovery of the axis ratio q, the potential parameters A and B using Monte 
Carlo simulations. The true values are q — 0.75, A — 1 and B — 2. The cases with 5% and 
15% errors on the flux ratios are shown. Notice that q and A are recovered well, but the 
larger (~ 15%) errors have a deleterious effect on recovery of B. The histograms show the 
distributions obtained with 200 realisations. 
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5.1. Tests for Elliptic Power-Law Lenses 

From an observational perspective, only the flux ratios of the images and not their 
absolute magnifications are measured. The flux ratios are defined as fi = So, they 

are normalised to the flux of the reference image /xi, which most observers generally take as 
the brightest image. For the quadruple lenses, there are three flux ratios /2, /s and /4, with 
/i = 1 by definition. So, the moments of the flux ratios are: 

4 4 

i=l i=l 

In other words, Fq = Fq^^ = -Po,y is just the magnification invariant normalised by the flux of 
the reference image /ii. Similarly, the higher moments of the flux ratios are given by dividing 
the configuration invariants by /ii. 

The source position is just 

e= (1 + 71)^ + 72^, r^ = 72^ + (l-7i)§^. (80) 

Substituting these into the first moments (33) and the reciprocal first moments (41), we 
obtain 

F-l,xFl^X ^ F-l^yFl^y 

Provided the lensing galaxy can be identified and the position angle of its major axis can 
be measured, then all quantities in eqn (81) are directly available from the data. So, these 
equations provide two simple, independent checks that are easy to apply. If they are satisifed, 
then the lens may be well-represented by an elliptic power-law potential with shear. If the 
major axis of the lensing galaxy is not known, then one of equations can be used to estimate 
the position angle and the second used to check that the elliptic power-law potential is a 
reasonable model. 

The main difficulty in application of the lensing invariants is the uncertainty in the 
observable quantities. The relative positions of the lensing galaxy and the images may be 
accurate to better than 1%. In propitious circumstances, the relative fluxes may be good to 
perhaps 5%. Unhappily, microlensing may render these fluxes much more uncertain, with 
the error perhaps even as high as 50% in the optical (Schechter 2000). Microlensing is of 
course valuable, as it sets constraints on the mass fraction in compact objects in the haloes 
of the lensing galaxies (e.g., Witt, Mao & Schechter 1995). However, for all researchers 
modelling the smooth, large-scale distribution of mass in the lensing galaxy, microlensing is 
a serious encumbrance. It is sometimes claimed that microlensing is primarily a difficulty 
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in the optical wavebands and may be circumvented using radio fluxes. However, there have 
been very recent identiflcations of microlensing events in the radio (Koopmans & de Bruyn 
2000), and so even the radio fluxes may be more uncertain than is widely believed. As 
benchmark values for our calculations, we take 5% and 15% errors. In some cases, the radio 
fluxes may already be accurate to 5%; in other cases, this accuracy will be achievable in the 
near future, once the effects of microlensing are better understood and quantifled. So, 15% 
is a conservative value for the errors in present-day radio flux measurements. 

Figure 3 shows histograms of the distribution of the ratios (81) in 400 cases, allowing 
for 5% and 15% errors in the flux measurements. In both cases, the distributions are peaked 
around unity, but there is substantial scatter. Even for datasets with flux errors of 5%, the 
test should be used in the form 

Ezl^Ihl^l±2, Ezl^EhlL^l±2. (82) 

for practical application. If the flux errors are as high as 15%, the tests become much less 
reliable. If a lens system passes this test (as for example the Einstein Cross does), then it is 
worthwhile to proceed to estimation of the lens parameters. 

5.2. Monte Ccirlo Simulations 

There are three galaxy model parameters, namely A, q and j3. In the general case, there 
are two unknown components of shear 71 and 72. Finally, there are three more unknowns, 
namely the source offsets ^ and t] and the position angle of the lensing galaxy. It is 
possible to solve formally for the model parameters, given the magniflcation invariant and 
the flrst, second and third moments. However, this is not a useful way too proceed as small 





-Aa 




Radio fluxes 


Radio flux 




(in ") 


(in ") 


(in 11 Jy) 


ratios 


A 


-0.09 


-0.94 


65.5 


1.00 


B 


0.58 


0.74 


64.2 


0.98 


C 


-0.72 


0.27 


26.5 


0.40 


D 


0.76 


-0.42 


59.4 


0.91 



Table 1: Observational data on the Einstein Cross. The positions of the images are relative 
to the centre of the lensing galaxy and taken from Crane et al. (1991). The radio fluxes are 
provided by Falco, as reported by Keeton & Kochanek (1996). 
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Parameter 


Best Fit Value 


A 


0.762 


Q 


0.885 




0.116 


7i 


0.117 


e 


67.25° 




-0.1062 


V 


-0.0266 



Table 2: Estimates of the model parameters for the lensing galaxy of the Einstein Cross. 



-Aa AS Flux 

(in ") (in ") ratios 

A -0.10 -0.92 1.00 

B 0.58 0.76 0.90 

C -0.72 0.26 0.40 

D 0.77 -0.40 0.83 



Table 3: Model fit to the Einstein Cross using the radio data. This can be compared with 
the data given in Table 1. 
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errors in the observables amplify into larger errors in the first, second and third moments. 
This is exacerbated by the plus and minus pattern imposed by the parities. The fitting of 
lens models to data normally proceeds by the minimisation of the differences between the 
observed and true positions (and flux ratios, if desired) using a downhill simplex routine 
to search parameter space (e.g., Kochanek 1991). The lensing invariants can be fruitfully 
exploited to reduce the dimensionality of the parameter space in which the search proceeds. 

As an example, let us consider a simple case. Suppose that we assume that the lensing 
galaxy has negligible external shear. First, it is straightfoward to find the source offsets 

C = -B > n = -B ■ (83) 

Prom the reciprocal second moments, it is then easy to see that the flattening of the potential 
is 



Prom the second moments, we obtain 



2B _ -^2,xi^o - Fl^ 
? - —ET-B (85j 



F2,yFo — -F\ 



y 



from which B can be estimated. There remain an unknown position angle of the lensing 
galaxy and an unknown potential normalisation, which can be obtained by minimizing 



{A, Q)= x^il - u2 + y2q-2)i-m )\ y^^^ ~ ix' + 'X-'y-f'/' ^] ^^^^ 

images \ i Hi'i I \ i yiH J 



Figure 3 shows the results of Monte Carlo simulations using fake data. Taking A = 1,5 = 2 
and q = 0.75, the lens equation is flrst solved to flnd the true image positions and fluxes 
for random source offsets. A random rotation is performed to take the image positions from 
the intinsic coordinate system of the lensing galaxy to the observer's coordinate system on 
the plane of the sky. The fluxes are adusted by random errors of 5% and 15% to give the 
simulated data. This is fed into the algorithm and the random position angle © and the 
parameters of the lensing galaxy A, B and q are recovered. Figure 4 shows histograms of 200 
such trials. We see that if the relative fluxes arc accurate to 5% or even 15%, then there is 
good recovery of the model parameters q and A. The parameter controlling the fall-off the 
potential B seems more vulnerable to errors, as there are tails to the distribution for the 
case of 15% errors. Tests also show that essentially no useful information can be extracted 
if the fluxes are only accurate to 50%. 
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5.3. The Einstein Cross 

Table 1 presents the available astrometric data on the Einstein Cross quadruple system 
from Crane et al. (1991). The astrometry is provided in terms of the relative right ascension 
and declination (in arcseconds) with respect to the center of the lensing galaxy. The radio 
data on the fluxes of the four images A,B,C and D is given, rather than the optical data 
which is known to be strongly effected by microlensing (e.g., Ostensen et al. 1994). The 
bulge of the lensing galaxy has a measured position angle of 68° on the sky (Rix et al. 1992). 
If the misalignment of the source is not too great, it is natural to expect the images of 
negative parity are located roughly parallel to the major axis, those of positive parity are 
roughly perpendicular to the major axis. Accordingly, we reckon that the negative parity 
images are C and D, whilst the positive parity ones are A and B. 

To fit an elliptical power-law model in the presence of external on-axis shear is straight- 
forward. We use the relations 



to minimise the 

X'(7l,AAe)= J] [e-X,(l+7i- ^ ^2g-2)l-/3/2 )] + [^7-^^(1-71- ^ ^2g-2)l-/3/2 )] 



images 



(88) 

over 7i, ^4, 0. The results for the model parameters are given in Table 2. First, let us note 
that the recovered position angle agrees well with the observed position angle of the major 
axis of the lensing galaxy. Second, only mild on-axis shear is required. The flattening of the 
lensing galaxy is significant, but not unreasonable. Third, the lensing galaxy is not close to 
isothermal, and there are no good fits with (3 = 1. As all the images are at almost the same 
projected distance, the density fall-off, and hence the parameter /3, is not well-constrained. 
The total mass enclosed within the images is more robust and comparable to the values 
found by previous investigators (e.g., Rix et al. 1992). Table 3 shows how well the best 
fit reproduces the data. The positions and the relative brightnesses of the four images are 
reasonably well reproduced. There are mild discrepancies in the magnifications of images B 
and D, which are both somewhat fainter in our model. The computational advantage of the 
lensing invariants (87) is that the minimisation can proceed in fewer dimensions, with 
obvious gains in accuracy and speed. 

Note that Witt & Mao (2000) give a cautionary example of a model with a flat rotation 
curve and with a projected density distribution stratified on ellipses. Although superficially 
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similar to the elliptic power-law models, the lensing invariants are quite different. For exam- 
ple the sum of the signed magnifications is ~ 2.8 rather than 2. The magnifications depend 
on the second derivatives of the lensing potential and so are more susceptible to the details of 
the modelling than the image positions. Real lenses certainly do exhibit diverse behaviour. 
Some - for example, B1422+231 where the lens is both a galaxy and a cluster or B1608+656 
where the lens is a pair of interacting galaxies - are not likely to be well described by any 
simple models. Nonetheless, much of the optical depth is in isolated giant eUiptical galaxies 
or in spirals dominated by dark haloes. It is reasonable to model such lenses with simple 
potentials and the algorithms we have developed have a real role to play. 



6. Conclusions 

This paper has presented a powerful new method for obtaining lensing invariants. Such 
invariants are calculable from the observables and they remain unchanged as the source 
of radiation is moved within the central caustic. The magnification invariant is the sum 
of the signed magnifications of the images, whereas the configuration invariants are sums 
over powers (either positive or negative) of the image positions multiplied by the signed 
magnifications. The idea behind the method is to use Cauchy's theorem and the residue 
calculus to recast sums over the images as contour integrals. 

We have illustrated the method by application to a simple family of galaxy models, 
the power-law models in which the gravitational potential is scale-free and stratified on 
similar concentric ellipsoids. They may be embedded in an external shear field of arbitrary 
orientation. The convergence k of the models falls off with distance like a power-law with 
index —2/B. The lensing invariants are exact for the models with B = 1 (the point mass 
case), B = 2 (an isothermal galaxy with projected potential ip oc r) and i? = 3 (a galaxy 
with ip oc r^/^). Some of the lensing invariants for the power-law models were previously 
calculated by Witt & Mao (2000) for the case of on-axis shear. Our results complement 
theirs by extending the invariants into the regime of arbitrary shear, by giving general and 
more extensive formulae, as well as by providing some entirely new invariants (the reciprocal 
moments) . 

In their discovery papers, Witt & Mao (1995, 2000) used direct elimination methods to 
find the first examples of lensing invariants. We believe that our contour integral method of- 
fers considerable advantages over this. All that is needed is the identification of the poles and 
the branch points of the lensing equation and the calculation of the residues using Cauchy's 
theorem. Our method is simpler than the mult i- dimensional residue calculus proposed by 
Dalai & Rabin (2000). Nevertheless, the idea motivating Dalai & Rabin's work - namely. 
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that sums over images may be replaced by contour integrals - is the same as that behind 
our paper. Only the familiar single-variable calculus is required in our formulation of the 
problem, but it has allowed us to go beyond the polynomial equations and meromorphic 
functions studied by Dalai & Rabin, and to handle branch points. In turn, this led us to 
our discovery of full series expansions for the lensing invariants of the four bright sources, 
unsullied by spurious roots of the imaging equation. 

We have shown how the lensing invariants may be exploited in modelling by using Monte 
Carlo simulations to demonstrate the feasibility of the recovery of parameters, such as the 
flattening and the normalisation of the potential well. Most lens fltting is conventionally 
done by minimisation. The lensing invariants play an important role by reducing the 
dimensionality of the space in which the global minimum is sought, with consequent gains 
in accuracy and speed. The power-law models with shear are sensible ones to use when the 
lens is a reasonably isolated spiral or elliptical galaxy. To illustrate this, we have used the 
lensing invariants to provide a good fit to the astrometry and photometry of the Einstein 
Cross (G2237+0305). 

Of course, caution is needed in application of the lensing invariants identified in this 
paper, as they only apply to a particular family of galaxy models. The pressing problem is 
to calculate the lensing invariants for more general potentials, such as galaxy models with 
cores or density distributions stratiflcd on ellipses (e.g. Kassiola & Kovner 1993, Barkana 
1998). This will lead to a greater understanding of how the invariants can constrain the mass 
distribution for real multiply-lenscd systems. We believe that our contour integral method 
can be adapted to these cases and are presently pursuing further investigations. 

We are indebted to Hans Witt for useful conversations and much encouragement. NWE 
thanks the Royal Society for financial support. CH thanks the sub-department of Theoretical 
Physics, University of Oxford for hospitality during a stimulating six-month sabbatical visit, 
and the Florida State University for awarding that sabbatical. His work is supported in part 
by NSF through grant DMS-9704615. 
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A. Images and Caustics 

When the source is sufficiently close to the center of the lens and \C\ small, then there are 
four images if S < 2 and five images ii B > 2. To see this, we split K into two components 
K — Ki + K2, as defined in eqn (57). We recall that, provided that 7i + 7I < 1, both ti and 
t2 (defined in eqn (34)) are positive with ti > t2 and ^2 < 1- As Figure 5 shows, the graph 
of —Ki, with its two zero minima at t — ti and t — t2, has two nearby pairs of intersections 
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with the positive K2, which is small when |C| is small. There is a fifth image for large positive 
t when B > 2 because, however small \(\, the asymptotic t^'''^|CP growth of K2 eventually 
overtakes the growth of —Ki. This fifth root is given by the n — case of equation (50). 

The number of images changes when pairs of roots of the imaging equation merge and 
become double. Both K and Kt vanish at double roots, and then 

1 dK2 1 dKi B 2 2 C C . .a.n 



K2 dt dt t P-\Q\ P + \Q\ PC-QC PC-QC 

Because each denominator is linear in t, equation (Al) gives a quartic in t, degenerating 
to a cubic in the special case of i? = 2. Moreover, the equation is independent of 
and depends only on the angular argument of the complex source position ( = |C|e"^. 
Hence the source positions which give double images are found by looking for positive roots 
of equation (Al) along each specific angular direction. We demonstrate below that, with 
our restrictions on the lens parameters, eqn (Al) has one real root in [^21^1] for all B, one 
real root in [ti,oo) for B > 2, and no others. Once t is found the imaging equation gives 
|(|2 = [p2 _ |g|2]2/[^i?(p _ Qe-2'<^)(P - Qc'^''^)]. The root in [t2,h] corresponds to the 
one point on the tangential caustic in the source direction 0, and the root in [ti,oo) for 
B > 2 corresponds to the one point on the radial caustic in that direction. These roots 
must generally be found numerically, but there are four special directions for which there are 
analytical formulas. These are the directions for which Qe"^"^ is real. When e^"^ = Q/\Q\: 
then t — t2 where P — \Q\ is a triple root of the quartic, and gives a cusp at |C| = 2|(5|/tf''^ 
on the tangential caustic. The remaining root is t — Bti/{B — 2), which for B > 2 gives 
a point Id = (2/S)[(5 - 2)/5ti](-^-2)/2 on the radial caustic. When e^'"^ = -Q/\Q\, then 
t = ti where P = —\Q\ is a triple root of the quartic, and t = Bt2/{B — 2) is the fourth 

B /2 

root. The former gives a cusp on the tangential caustic at which \(\ = 2\Q\/t-^ for B < 2, 
and also for i? > 2 if ti is the smaller of the roots. In the latter case, the fourth root gives a 
point Id = {2/B)[{B - 2)/Bt2Y^-^^/^ on the radial caustic. However, if Bt2/{B - 2) < h, 
then the latter point is on the tangential caustic, while |d = 2|(5|/ti on the radial one. 
Depending on the parameters, either caustic may lie wholly inside the other, or they may 
intersect as in the case shown in Figure 6. 

In the general case, QC/C = IQ|e'^ for some angle Q which is not an integer multiple of 
TT. To find the roots of equation (Al) for —1 < cos^ < 1, we rewrite that equation as 

B_ 2(P-|g|cos^) 2 2 _ -2p3 + 6p2|g|cos^-6P|g|2 + 2|g|3cos^ 

7 ~ p2+|g|2-2P|g| cos^ ^ t^2 ~ (p2-|g|2)(p2+|g|2_2P|g|cos^) ■ 

(A2) 

Graphs of its two sides as functions of t are shown in Figure 7. The right hand side has 
vertical asymptotes at t = ^2 and t — and a zero in between. This is the only zero of 
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Fig. 5. — Graphs of K2{t) (dashed curve) and —Ki{t) (full curve) for the case B — 3, q — 0.8, 
no shear and C = 0.15 exp(i7r/3). The two curves intersect a fifth time at a much larger value 
of t. The intersections correspond to the five images. 




Fig. 6. — The two caustics for the lens of Figure 5. Note that the cusps of the tangential 
caustic are "naked" and protrude beyond the radial caustic (see e.g., Evans & Wilkinson 
(1998)). 
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Fig. 7. — B/t (dashed curve) and the right hand side of equation (A2) (full curve) for the 
case S = 3, g = 0.8, and no shear. The caustics are shown in Figure 6. The angle 9 = ti/?>. 
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the numerator because its derivative with respect to t, which is the negative of that with 
respect to P, is 6[P — \Q\ cos^]^ + 6|Qpsin^^ and hence always positive. The two sides of 
equation (A2) have different signs and hence no roots in (0,t2))- They have at least one 
root in {t2,ti), and another in (ti, oo) for B > 2 because the right hand side, which decays 
as 2/t is ultimately the smaller for large t. To prove that these are the only possibilities, 
we rewrite equation (A2) as a quartic in two different ways, and apply Descartes' rule of 
signs (Henrici 1974). We also reduce the number of parameters to three by working with the 
variable a — {t — ti)/\Q\ — —1 — P/\Q\. The first form of the quartic is 

{B-2)a^ + 2[{B -2) {2 + cos 9) - {po + cos 9)]a^ 

+ 2(1 + cos 6) [3{B -po- 3)a^ + 2{B - 3po ~ 5)a - 4(1 + po)] = 0, (A3) 

where po = Po/\Q\ > 1. Let c„ denotes the coefficient of c". Each coefficient is manifestly 
either zero or negative for B < 2, and hence there can then be no roots with a > i.e., t > ti. 
For i? > 2, C4 > 0, Co < 0, and C2 - Ci = 2(1 + cos9){B + 1 + 3po) > and so C2 > Ci. We 
now show that, regardless of the sign of C2, there can be only one sign change from positive 
to negative in the sequence of coefficients (04, C3, C2, ci, cq). If C2 > 0, then B > po + 3 and 
C3 > 4 + 2po(l + cos^) > 0, and the single sign change occurs somewhere between C2 and 
Cq. If C2 < 0, there is a single sign change somewhere between C4 and C2, depending on the 
sign of C3. By Descartes' rule of signs, the single sign change in the sequence of coefficients 
means that the quartic has exactly one root in > 1. 

To show that there is one and only one root of the quartic in (^1,^2), we use a bilinear 
transformation to the variable r = {t — t2)/{ti — t) — {\Q\ — P)/{\Q\ + P) which maps the 
interval ^2 < i < to the whole of r > 0. Equation (A2) remains a quartic, and becomes 



(1 + cos ^) [(1 + poy + {B+po- ly] - (1 - cos ^) [(1 + Po - B)t + {po - 1)] = 0. (A4) 



Now C4 > 0, C3 > 0, C2 = 0, and Cq < 0, when c„ now denotes the coefficient of r"^. There 
is a single sign change from positive to negative in the coefficient sequence regardless of the 
sign of Ci, and hence, by Descartes' rule of signs, always one root for t in (^2,^1)- This root 
gives a point on the tangential caustic. 



Witt & Mao (2000) give the third order moments in the case of on-axis shear only. The 
third order moments, derived from the residues of simple poles at t — 0, and double poles at 
t — ti and t — t2 are 



B. Third Order Moments 




(Bl) 
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E 



q^j , , i3[(l-7i)e-72r/]^[(l + 7i)^-72e] 

4' (1-7? -72)' 



^[(l-7i)e-72^][(l + 7i)^-72e]' 

(1 - 7? - iir 



(B2) 
(B3) 

(l-7?-72Y ■ ^^'^ 
The Ri, R2, /i, and I2 terms are all linear in the source coordinates. They are the real and 
imaginary terms of the m = 3, n = and m = 2, n = 1 configuration moments respectively, 
and are 



3 q\oT , i^[(l + 7l)^-72e]- 

2^ mi = 7(3^2-/1) + — 7^3-J3-|r^ 



R2 = 



''1 



2 + 



1 - 



{B + 1)\Q\ 



1 + 



{B + l)\Q\ 

t2 



{B + 1)\Q\ 
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''2 



2 - 



iB + l)\Q\ 

t2 



X 



[(Re(g))2e + 272?'Re(g)77 - 7!?'^] 



iB + l)\Q\i 
1 



fB+2 



fB+2 



+ 



' ^B+2 
L2 t-i 



+ 



\Q\ \tt 



^B+l ^B+1 



+ 



1 



j.B+1 



2 + 



1 - 



1 + 



[Re(g)e + 72g'r;]l>, 

(i?+i)iQr 



[Re(g)7y + 72^] 



{B + 1)\Q\ 



+ 



^B+l 



{B + 1)\Q\ 



[{Re{Q)frj - 2^2MQ)^ - iWn] 



q^B 



+ 



[B + l) 



\Q\ 

B + i)\Q\v 
1 



^B+2 



fB+2 



+ 



.B+2 ' .B+2 
12 f-l 



+ 



\Q\ Vtf+^ t^2^' 



[Re(g)r; - 72^] 



(B5) 



C. Formulae for Hyper geometric Functions 

The hypergeometric functions that occur in this paper are all of the form 2-Fi(M, jB; M+ 
N]l — z), where M and N are integers and both are positive. The following formula expresses 
them as two finite sums, one with N terms and the other with M terms. It is obtained by 
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applying first Abramowitz and Stegun (1965, hereafter AS) equation (15.3.9) and then their 
equation (15.3.5). The result is 

(M+Ar-l)!^ ' ' ^ iN-^)^t^]tSl:t^ijB-s) 

'' A:=0 11 



It is significant that the products of terms in B in the denominators are all subparts of a 
longer product of terms in B which multiplies the hypcrgcomctric functions in the definition 
(74) of the coefficients /j,£,Ar which appear in our series for the lensing invariants. Hence, 
these coefficients never have poles in B and are all analytic functions of B. For example, the 
two N — integrals needed for the evaluation of the extra term in equation (66) are: 



h,2 



Tern 



-(S-2)(S-3)+2(S-l)(S-3) ( -1 l-(S-l)(S-2) ( ^ 1 +2 



f2,B-3 
''V'2 



[B-l){B-2)-2{B-l){B-3) (^f^ +{B-2){B-3) (^fj -2 (^fj 



2 \ 3-B 

i6|g|3 

The similarity of these two expressions is an instance of the general symmetry relation: 

Ij,e,N{tl: t2) = Ij,2j-e,N{t2: h). (C2) 

Hence, only two distinct / terms are needed for the first order corrections for the first 
moments 

E ^'^p^^^ - ^^Iy^^^^^+kI^^^ (C3) 

4 images ^ ll I2) 



4 images 



(l-7?-7i)^ 



where 

= b^Ii,o,i + abIi,i,i + abIi^2,i + a^h,3,i, 

A-io ^ cg[-36/i,o,i + (6-2a)/i,i,i + (26-a)/i,2,i + 3a/i,3,i], 

= g'[(a6+2c2)Ji,o,i + (a'-2c2)/i,i,i + (62-2c2)/i,2,i+(a6+2c2)/i,3,i], 

^03 = cq^[-ah,o,i + ah,i,i-bIi,2,i + bIi,3,i\, 

Kq3 = eg [-6/1,0,1 + 6/1,1,1-0/1,2,1 + a/1,3,1], 

i^oj = g2[(a^+2c^)/i,o,i + (6'-2c^)/i,i,i + (a^-2c2)/i,2,i + (a6+2c^)/i,3,i], 

KI2 = cg3[-3a/i,o,i + (a-26)/i,i,i + (2a-6)/i,2,i + 36/i,3,i], 

Kl^ = g'[a2/i,o,i+a6/i,i,i+a6/i,2,i+6'/i,3,i], (C5) 
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and the constants a, 6, and c are 



a — 



Bq' 



Reg 



6 = 



Bq' 



Reg 



c — 



2|g| 



(C6) 



Finally, let us caution against using eq. (CI) for numerical calculations for all but small 
values of M and iV, and for z close to 1, because it can become a difference between two 
nearly equal quantities. The hypergeometric series, which converges rapidly for z close to 1, 
should be used instead. 



D. Convergence of Series and Analytic Continuation 

To prove that the series of equation (65) converges for sufficiently small we need a 
bound on the coefficients Ij^2m.- The cases B <2 and B>2 require different treatments. For 
B <2, we use the integral (63), replace the powers of (w + ti) in the denominator with the 
smaller {w+t2), and derive the bound 

1 r w^'^'dw AB-2)j-2 r v^^'-'dv 



The simpler integral can be evaluated in terms of Gamma functions (AS, eq. 6.2.1). As this 
bound is independent of m, we can take it out of the inner summation over m when we use it 
in the infinite series in equation (65). The inner sum can then be carried out to give simply 
(|Ap + |z/p)-' = ICp-'. The series for the magnification is therefore majorized by the series 

Bq' ^ r[ji?]r[2 + (2-i?)j] , 

(2^Tl)! ' ■ ^ ^ 

Applying the ratio test and using Stirling's formula (AS, eq. 6.1.39) to approximate the 
Gamma functions, we find that the series (65) converges if 

Id < 2 j B-^/'. (D3) 

This condition is met uniformly for all B in [1, 2] when 

ICI < v^, (D4) 
and then the infinite series (65) for M.{B) is absolutely and uniformly convergent. 

The integral (63) is not valid for B>2, and we must then use equation (64) and find 
a bound for the hypergeometric function. We can use the hypergeometric series (AS, eq. 
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15.1.1) for that function because its fourth argument hes in (0, 1). As all the other arguments 
are positive, all the terms of the infinite series are positive. Moreover, each of them is less 
than the corresponding term of the infinite series for 2-^1 (2j + 2, ji?; 2j + 2; 1 — ^2/^1) because 
a. < 2j. This majorizing series is geometric (AS, eq. 15.1.8) and sums to {t2/ti)~^^. The 
terms (jB — s) in the product are now all positive for large enough j because s < 2j + l, so 
we ignore the s terms and obtain, for large j, the bound 



< 



jjB)^^^ ^2mr-2j-l^jB-2mr-l 



ti 



(2j + l)! 2 

The series for the magnification is now majorized by the series 



tit2 



E 



{JBf^^ 
(27 + 1)! ^ 

^ ' ^ m=0 



2m 



2m-2j 



The inner sum is 



^2 



(D5) 



(D6) 



(D7) 



and its maximum possible value is ICPVV because |Ap + |z/p = Hence, the infinite 

series (65) for the magnification invariant is majorized for S > 2 by the series 



5g2 

tit2 



^{jBy^^H^ f\c\ 



(2j + l)! 



f2 



(D8) 



The ratio test show that this series is convergent, and that the infinite series for the magni- 
fication is therefore absolutely convergent, if 



ICI< 



2U 



eBt\ 



B/2 ■ 



(D9) 



Thus, however large B is, the series for the magnification invariant converges for sufficiently 
small offset. Actually, if ti, the larger root of the equation = is less than unity (as 
it is if the on-axis shear 71 is negative) and 



then 



-71 > 



Btf'^ < 



l-g2 



In 



(DIO) 



(Dll) 



for all B, and the magnification series (65) has a finite radius of convergence that is inde- 
pendent of B. When ti — 1 {as it is when there is no shear), then the bound (D9) decreases 
as 1/B. 
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Regardless, if we need the magnification invariant M.{B) for a model for which B — 
-Smax > 2, then we can find a S, independent of B, such that the series for Ai{B) converges 
absolutely for |(^| < 5 for 1 < -B < -Bmax- That is, our series is uniformly convergent as a 
function of B. The explicit expressions that we have for each term in the series (65), from 
equation (64) and equation (CI) of Appendix C, shows that each is composed of polynomials 
and exponential functions of B. Hence, each term of the uniformly convergent series is an 
analytic function of B. It follows, by a well-known theorem of complex analysis (e.g Levinson 
& Redheffer 1970), that the sum M{B) of the series is also an analytic function of B for 
Id < S. It can therefore be continued beyond the range 1 < S < 2 for which it was originally 
derived by the analysis of §4.1 to apply to models with larger B. 

With only minor changes, the series for the configuration invariants can be handled 
in a similar way with a similar outcome; the infinite series (73) for Sm,n{B) is absolutely 
convergent when the same inequahties (D4) and (D9) as for M.{B) hold. Hence, the series 
for the configuration invariants of the four bright images also remain valid for B > 2 for 
sufficiently small offset \(\. 

Formulas for one of the four bright images become singular at the radial caustic where it 
merges with the weak fifth image. Therefore we can not expect our series to converge there, 
and inequality (D9) does not hold at the known analytical \(\ values given in Appendix A 
for special points on the radial caustic. However, the tangential caustic need not limit the 
convergence, even though series lose physical significance there, because the singularities of 
the two merging images can cancel. One can verify that the B < 2 inequality (D3) may 
indeed hold at the cusps of the tangential caustic. 



